PULSATILE VISCOUS FLOWS IN ELLIPTICAL VESSELS AND ANNULI: 
SOLUTION TO THE INVERSE PROBLEM, WITH APPLICATION TO BLOOD AND 
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Abstract. We consider the fully-developed flow of an incompressible Newtonian fluid in a cylindrical vessel with 
elliptical cross-section (both an ellipse and the annulus between two confocal ellipses). In particular, we address an inverse 
problem, namely to compute the velocity field associated with a given, time-periodic flow rate. This is motivated by the 
fact that flow rate is the main physical quantity which can be actually measured in many practical situations. We propose 
a novel numerical strategy, which is nonetheless grounded on several analytical relations. The proposed method leads to 
the solution of some simple ordinary differential systems. It holds promise to be more amenable to implementation than 
previous approaches, which are substantially based on the challenging computation of Mathieu functions. Some numerical 
results are reported, based on measured data for human blood flow in the internal carotid artery, and cerebrospinal fluid 
(CSF) flow in the upper cervical region of the human spine. As expected, computational efficiency is the main asset of our 
solution: a speed-up factor over 10 3 was obtained, compared to more elaborate numerical approaches. 

The main goal of the present study is to provide an improved source of initial/boundary data for more ambitious 
numerical approaches, as well as a benchmark solution for pulsatile flows in elliptical sections with given flow rate. The 
proposed method can be effectively applied to bio-fluid dynamics investigations (possibly addressing key aspects of relevant 
diseases), to biomedical applications (including targeted drug delivery and energy harvesting for implantable devices), up 
to longer-term medical microrobotics applications. 

1. Introduction. Pulsatile flows, here intended as flows driven by a time-periodic force (generally 
the pressure gradient), are encountered in many contexts of practical interest. Amongst them, a remark- 
able case is that of heart-beat driven, human physiological flows, including blood circulation [50] and 
(less directly) cerebrospinal fluid (CSF) flow (3TJ[57j- Besides bio-fluid dynamics, time-periodic flows are 
widely studied with regard to e.g. chemical-physics applications, like species separation [30] . as well as 
for mass and heat transfer problems [331 158) . Furthermore, many studies can be found in the literature 
which address pulsatile flows in a more general (yet still application-driven) fluidic context, like peri- 
staltic pumping [46] . to name but one. However, in many cases of practical interest the pressure gradient 
is unknown or hardly measurable, and the only available datum is the time-dependent flux (hereafter, 
flux and flow rate are understood as synonyms). This is the most common measurement, for instance, 
in clinical settings regarding blood and CSF flow. 

As for blood flow, it is difficult for many portions of the vasculature to keep an ideal, circular 
cross-section, due to the presence of surrounding organs (possibly activated by relevant muscles, or 
simply displaced by postural changes). Hence, while keeping some degree of abstractness, it is more 
appropriate to consider an elliptical cross-section to obtain a more accurate description of the flow (see 
e.g. [H|). As for CSF flow, the spinal sub-arachnoid space can be well approximated by an elliptical 
annulus [31 j : CSF dynamics in such a domain is affected by pulsatility (mainly towards the upper 
cervical region) and plays a major role in the (still poorly understood) pathophysiology of many high- 
impact disease, like syringomyelia and Chiari malformation |T5l [38] . Of course, in order to tackle realistic 
3D flow conditions, for instance associated with patient-specific geometries, a fully numerical approach 
is mandatory (see e.g. [T51 [2H1 ESI 123 for the CSF case). This implies the exploitation of rather time- 
consuming approaches, even when adopting a classical Newtonian behavior for the working fluid (as 
surely appropriate for CSF [37], as well as for blood in large vessels [TT1 [48] ) ■ Nevertheless, it is possible 
to keep some degree of physical representativeness also by adopting the hypothesis of a fully-developed 
flow (see e.g. [24] ). despite its inherent oversimplifications from a physical viewpoint. Such a position 
has been systematically adopted in many studies in order to obtain analytical flow solutions (see below), 
otherwise hardly achievable. However, even by considering a fully-developed flow, it is necessary to solve 
a non-standard, inverse problem in order to evaluate the velocity field (and the pressure gradient) by 
starting from a given flux across any cross-section of a vessel. 

With reference to the above observations, we point out to expressly confine our attention to a 
special class of flow solutions, namely fully-developed flows of an incompressible Newtonian fluid in a 
straight cylinder with elliptical cross-section (either simply connected or not, see below). In particular, 
these assumptions permit to address a simplified linear problem, for which superposition principle holds. 
Nevertheless, the adopted positions aim at obtaining a benchmark solution for the inverse problem 
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described below: this is the main target of the present work. More precisely, by virtue of the fully- 
developed flow hypothesis, we provide a numerical strategy which is substantially grounded on analytical 
relations, while being more amenable to implementation than previous ones proposed in literature. Hence, 
while laying no claims of generality (not even deliberately addressing realistic three-dimensional flows, 
e.g. within patient specific geometries and possibly involving complex rheological effects), the method 
we propose is credited to hold some value for obtaining at least an approximation to real flow conditions, 
while containing the computational cost, in the spirit of benchmark flow solutions. Furthermore, the 
obtained solution strategy can be used as an improved source of boundary data for more ambitious 
numerical approaches based on realistic data, up to serving as a debugging tool, ft is worth stressing 
once again that we address the inverse problem, for which no straightforward approaches have been well 
established so far. 

To better define the problem, we consider the motion of an incompressible (the constant density is 
here normalized to 1), Newtonian fluid in a semi-infinite straight pipe P = E x R + c R 3 . In the sequel 
E C R 2 will be cither an ellipse or an elliptical annulus (even if results for circular cross-sections are also 
recalled, and some abstract results can be obtained for any smooth and bounded cross-section). The 
governing equations are therefore the incompressible Navier-Stokes equations which, in a reference frame 
with z directed along the pipe axis (and x := (x\,X2) belonging to an orthogonal plane), read 

d t u + (u ■ V) u - vAu + Vp = (x,z) e E x R+, t G R + , 
V-u = (x,z) E E x R+, t E R+, 

u = Q (i,z)effixR + ,ieR+ 

where u(t, x, z) and p(t, x, z) respectively denote velocity and pressure, and v represents kinematic viscos- 
ity. Furthermore, as anticipated, we deliberately neglect turbulent effects by looking for a special class of 
solutions, namely fully-developed solutions (also named Poiseuille-type solutions) of the following form: 

u(t, x, z) = (0, 0, w(t, x)) and p{t,x,z) = — X(t,x,z) + po(t), 

where Po{t) denotes an arbitrary function which disappears from the equations, since pressure only enters 
through its spatial gradient in the formulation of the problem. In addition, the following flux condition 
is assumed to be satisfied: 

/ w(t,x)dx = /(<), 
Je 

for some given scalar function f(t), only depending on the time variable. By standard arguments, the 
above Poiseuille-type ansatz implies that we may assume that the pressure has the form p(t, z) — — X(t) z. 
Moreover, the dependence of w on the space variables x\ and x-i allows us to consider a problem reduced 
to the cross-section E, with the given flux condition. This implies that we have to study the following 
problem: Find (w(t,x),X(t)) such that 

' d t w{t, x) - vA x w(t,x) = X(t), x £ E, te R + , 

(f.i) <w(t,x) = xedE, teR + , 

J E w(t,x)dx = f(t) teR+, 

where A x denotes the Laplacian with respect to the variables x% and X2- This can be considered as an 
inverse problem since, if the axial pressure gradient A is known, one can evaluate w (by solving a two- 
dimensional scalar heat equation) and thus the corresponding flux as well. Contrarily, in our problem the 
datum is the flux f(t) and one wants to find the corresponding velocity and pressure associated with the 
fully-developed solution (such a problem is linked to one of the nowadays classical Leray's problems [2CJ]). 
For the sake of completeness, we recall that existence and uniqueness of the fully-developed solution in 
pipes with very general cross-sections (and with explicit but not exact expression in terms of the flux) 
have been recently obtained in [21 [TIES], under reasonable regularity on the periodic f(t). 

The relevance of the cross-section shape on the solution process deserves some discussion at this 
stage. Indeed, when considering a circular cross-section, one can write an explicit solution in terms of 
special functions. In particular, the analytical solution of the direct problem was firstly obtained by 
Sexl |51j and Womersley |59j in terms of Bessel functions, while we recently derived in [6] an analytical 
solution to the inverse problem involving regularized confluent hyper-geometric functions. Such a fully 
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analytical approach is made possible by the particularly simple expression taken by the Laplace operator 
in cylindrical coordinates. In fact, when considering an elliptical cross-section, in spite of the apparent 
similarity with the circular case, the situation drastically changes: an explicit solution can only be 
obtained for the stationary problem, while it is necessary to resort to numerical computations already 
when addressing the direct time-dependent problem |3H E3 [S5J |5B] . In particular, current solutions 
in literature are based on involved tools like the Mathieu functions 0D]. Nevertheless, some numerical 
approaches which also tackle the inverse problem have been recently proposed [21] , still by invoking the 
Mathieu functions, also involving the determination of the Laplacian eigenvalues in an elliptical domain. 
Yet this is not a trivial task: there are several works emphasizing how instabilities can arise in the 
numerical computations when addressing such a problem by directly managing the Mathieu functions, 
especially when the ellipse eccentricity is small |26j . Indeed, computing the Mathieu functions has been 
recognized as the most critical issue in |24| ; major efforts have been spent in order to define a wise 
algorithm for their evaluation. 

In the light of these points, we propose an alternative numerical approach to the inverse problem of 
a pulsatile fully-developed flow in elliptical vessels. Our approach only indirectly addresses the Mathieu 
problem and therefore it does not suffer from the aforementioned limitations. Some numerical simulations 
are also reported, based on flow rates actually measured for human blood flow in the internal carotid 
artery, and CSF flow in the upper cervical region of the human spine. Corresponding results support the 
effective usability of our formulation (successfully coping with small eccentricities as well), by virtue of 
the dramatic reduction in computational time it allows, as compared to more elaborated finite clement 
approaches. This, in turn, broadens the application spectrum of our solution, as highlighted in the last 
section of the present work. 

Plan of the paper: In Section [2] we review relevant issues regarding stationary, fully-developed 
flows in elliptical cross-sections, either simply connected or not. Some results related to circular cross- 
sections are also recalled, for completeness. In Section [3] we propose a novel, computationally-efficient 
numerical strategy for solving the inverse time-dependent problem, by also exploiting an auxiliary, di- 
rect formulation. In particular, besides recalling recent results for the the circular cross-section, we 
discuss in detail the case of the simply connected elliptical section (being the annular case very similar). 
Furthermore, exemplificative numerical simulations are reported in Section |4j based on blood and CSF 
physiological flow data; computational cost of the proposed strategy is also compared to more elaborate, 
finite element numerical approaches. Concluding remarks are finally reported in Section [5j where the 
potential for effective application of our solution is highlighted, with explicit references to the biomedical 
and microrobotic research field. 

2. Stationary problem. Despite being classical, the stationary solution in the elliptical cross- 
section is re-derived in a way that is functional to subsequent treatment. Classical results for the circular 
cross-section are recalled as well, for completeness. Such a stationary problem is formulated as a Poisson 
problem with Dirichlet boundary conditions (besides the flux condition), namely: Find (w(x),X) such 
that 

— vA x w(x) — A, x G E, 

(2.1) \ w{x) = 0, x G BE, 

J E w{x)dx = /, 

where / € R is the given flux, while (—A) G R is the constant pressure gradient. 

2.1. Stationary flow in a circular cross-section. The history of this problem dates back to the 
work of Hagen and Poiseuille (see e.g. [2U] for details), who firstly derived the solution in a circular domain. 
Nonetheless, the interest for such a special solution (which can be easily derived thanks to complete 
separation of variables) is still considerable, as a benchmark flow also in applications to turbulence (see 
e.g. [H]), or for innovative biomedical applications (like magnetic particle targeting [57]). This is mainly 
due to the fact that Poiseuille flow represents one of the few examples of exact solution of fluid equations 
which are not space-periodic, yet with Dirichlet boundary conditions. 

2.1.1. Flow in the circle. Let E be a circle with radius R > and centered at the origin, namely: 
E = {(xi,X2) € R 2 : \x\ < i?}, where \x\ = \J x\ + x\. For such a case the solution to (2.1) is easily 
obtained, namely the well-known parabolic profile 



2.1.2. Flow in the circular annulus. Let E be a circular annulus delimited by radii R\ and 
i?2 > Ri, namely: E := {(£1,2:2) G R 2 : Ri < \x\ < R2} ■ In this case the solution of the Poisson 
problem (2.1 ) is 



(2.3) 



w(x) 
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Notice the limiting behavior of ( 2.3 1 for R2 = R and i?i — >• 



2.2. Stationary flow in an elliptical cross-section. We review some well-known results, which 
will be needed later on to handle the time-dependent case. Relevant solutions for this problem originally 
appeared in the pioneering works of Khamrui [32] and Verma ES1 EH] (besides earlier, preliminary 
contributions cited therein). Also in the present case it is possible to obtain an analytical solution, thanks 



to the particularly simple expression of the external force in the Poisson problem (2.1 ). Nonetheless, it is 



worth mentioning that numerical approximation of the aforementioned problem in an elliptical domain, 
yet with more general external forces, is the subject of ongoing research (see e.g. |35j). 

2.2.1. Flow in the ellipse. Let E = {(xi,x 2 ) 6 R 2 : xf/a 2 + x\ //3 2 < 1} denote an ellipse, where 
a = a cosh(6) and f3 = a sinh(fr) respectively denote the length of the major (xi-direction) and minor 
semi-axis (^-direction), while 2a represents the inter-focal distance. We preliminarily observe that the 
equality 



(2.4) 



cosh 



= b = sinh 



from which b 



implies log (a /a + y/a 2 — a 2 /a) — log (j3/a + \/ (3 2 + a 2 /aj , 
and a = \J a 2 — f3 2 . Then, to construct the solution we introduce the natural change of variables 



(2.5) 

with Jacobian 
(2.6) 



x\ = acosh(?7) cos(6 l ), 
X2 = asinh(?7) sin(#), 



J(r),9) = a 2 [sinh 2 (rj) + sin 2 (0)] 



a 

y 



-2iB 



■ cosh(2?7) 



By employing this change of variables we can reduce the problem in the rectangular domain 
(2.7) E' = { (77, 9) G R 2 : < 77 < b and < 6 < 2tt} , 

and it is well-known that the Laplace operator is separable with these coordinates. Furthermore, by 



applying the aforementioned change of variables, the Poisson equation (|2.1|) becomes 

v I d 2 u 




= A 



J(r),6) \drj 2 

u is 27r-periodic in 9 

u(b,9) = 

u(r/,9)J{r),9)dr)d9 = f 



9 G [0,2tt[, 



where u(rj,9) simply represents w(x) in the new variables (please also notice that J > for rj > 0). 
Moreover, we observe that an additional Neumann condition is needed at 77 = 0, as detailed below. Due 
to periodicity in 9, the solution can be written with the following complex Fourier series expansion with 
respect to 9 (hereafter "hat" ? denotes a complex Fourier coefficient): 



neZ 



where 
1 



1 

27- 



u(r],9) e 



By virtue of linearity, we temporarily assume A = 1 and we recast the problem at hand as follows (for 
the function u which is 27r-periodic in 9 and temporarily we do not assign the flux): 



■"( 



d 2 u 0"u^ 

u(b,9)=0 9e[0,2n[. 



l)=J(v,0) (v,0)eE', 



We then plug in the Poisson equation the Fourier series for both u and J(r), 9), as given in (|2.6[ ), and we 
observe that, in order to have a real-valued solution, we need to impose as usual m_„ = u n (hereafter 
"bar" 7 denotes complex conjugacy). In particular, ^S(u ) = (hereafter 3? and 9 respectively denote the 
real and imaginary part). Moreover, by using the symmetries of the solution (reflections over the ellipse 
axes of symmetry) it follows that u n — for n odd. Furthermore, in order to have a smooth solution 
also at r\ = 0, we need to impose d v u(0,9) = 0. As a result, by exploiting the above conditions, and by 
observing that J only possesses three non-zero modes (namely those associated with n = —2,0,2), we 
obtain the following identity: 



n=-2fl.2 



cosh(2?7) 



which corresponds to a system of three complex ordinary differential equations. In the above expression 
and in the sequel, "prime" ( • ') denotes derivative with respect to rj. Then, by equating the terms 
corresponding to the same exponential, and by recalling the conditions at rj = 0, b coming from the 
Dirichlet condition and from the symmetries of the solution, we have to solve the following boundary 
value problems: 



S±a(»7) - 4 « 



±2 



4i/ ' 



and < ~ uv " 1v 

u (6) =u' (Q) = 



cosh(2?7), 



u ±2 (b) = »(«4 2 (0)) = 9f(u± 2 (0)) = 0, 
Hence, by explicitly solving these uncoupled linear differential equations we get 



Uo(r]) = — — (cosh(2?7) — cosh(26)) . 



a 2 1 



a 2b+2r) 



16u 



1 



o4b 



Then, by summing-up the solutions after multiplication by the corresponding complex exponential, and 
after some algebraic manipulations, the solution to the problem with A = 1 turns out to be 

2 

u(r},9) = (cosh(2fe) - cos(26»)) (cosh(26) - cosh(27y)) sech(26), 
of 



being the corresponding flux 



u(r),9)J(ri,9)dr]d9 = a 2 



u(t),0) [sinlr (t?) + sin 2 (9)] d9 dr) 



Z2v 



sinh 2 (26) tanh(26). 



That said, the solution to the original problem (2.1 ) with given flux / is finally obtained as (Au(ry, 6*), A), 
i.e. by scaling the function u(j], 9) computed right above by the following factor: 

(2.8) 



A = / [sinh 2 (26) tanh(2t 



accounting for the actual flow rate. Finally, by mapping back to Cartesian coordinates, the sought 
solution reads 



(2.9) 



w(x 1 ,x 2 ) = 



2/ 

Tra/3 



1 



It is interesting to note how the solution (2.9| derives from Poiseuille solution (2.2 1 by a simple anisotropic 
scaling of the variables (as for obtaining the considered elliptical domain by starting from the unit circle). 
Thus, the elegant expression ( |2.9| ) (appearing in e.g. [55], yet missing in many related works) could have 
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been obtained through simpler derivations. Nevertheless, we deliberately decided to introduce many 
details above, since they serve as useful references for the solution of the corresponding time-dependent 
problem. On regard, it should be noticed how, despite separation of variables (in elliptical coordinates), 
ellipticity per se creates a new situation already when considering stationary conditions: there appear 
3 coupled modes in the solution (contrarily to the fully decoupled flow solution in the circle). Such a 
difference persists and it is somehow magnified in the time-dependent case, where more complex behaviors 
naturally appear (see Section [3]). 

2.2.2. Flow in the elliptical annulus. Given the inter-focal distance 2a, it is straightforward to 
apply ( 2.5 ) in order to also map two confocal ellipses into a rectangular region. More precisely, given semi- 
axes ai,f)i and a 2 , (3 2 , with f3 2 > /3i, the rectangle E' — {(77, 9) £ R 2 : b\ < 77 < b 2 and < 9 < 27r} 
is obtained, where b± and b 2 > b\ are derived from the respective semi-axes through (2.4). For this 



problem we have Dirichlet boundary conditions at both 77 = b\ and 77 = b 2 . Hence, by enforcing the 
symmetries u_„ = u n needed to have real-valued velocities, we end up with the following boundary value 
problems (for the auxiliary problem with A = 1): 



±2 



a 

4z/ 



u± 2 (bi) = u± 2 (b 2 ) = 0, 



and 



Sofa) 



2v 

{»0 (61) = Uo(&2 



cosh(27y) 
0, 



whose solutions, by direct integration, read: 

- 2 r (f? - bi) cosh(26 2 ) - (t? - b 2 ) cosh(26i) 



u± 2 {v) 



(b 2 ~ 60 

e 2(77-6 2 ) _ e -2(r,-b 2 ) _ e 2(»j-&i) + e -2(ij-6i) 



e 2(b 1 -b 2 ) _ e -2(bx-b 2 ) 



— cosh(2r7) 
1 



so that the following stationary solution for the auxiliary problem is obtained (cf. |55|): 
" 2 r (77 - h) cosh(26 2 ) - (tj - b 2 ) cosh(26i) 



u(r),6) 



sinh(2(r7 - 61)) - sinh(2(r/ - b 2 )) 
sinh(6 2 — bi) 



cos(26») - (cosh(2?7) + cos(20)) 



Finally, by reasoning as in Section 2.2.1 the solution to (2.1 ) with given flux / is obtained as (Xu(r], 9), A), 
with (cf. [24] ) 



A = / 



16f 



(sinh(46 2 ) - sinh(46i)) (cosh(26 2 ) - cosh(26i)) 2 cosh(2(6 2 - h)) - 1 



2(6a-6i) 



sinh(2(& 2 - 61)) 



Following |55j , it is worth remarking that the solution in the (simply connected) ellipse cannot be obtained 
as a limit case of the one in the annulus, due to the constraint of confocality. Indeed, by allowing b\ — s- 0, 
the inner boundary of the annulus (where the Dirichlet condition also applies) degenerates to the inter- 
focal segment. The resulting flow field is therefore that one associated with an ellipse with semi-axes a 2 
and j3 2 , yet also containing a plate having width equal to the inter- focal distance 2a. 

3. Time-dependent problem. We now address the time-periodic elliptical case, by proposing a 
novel numerical method for its solution (this is the main asset of the present study). More in detail, we 
firstly recall the solution of the inverse problem for the flow in the circle, so as to keep some degree of 
symmetry with Section [2j We then address the elliptical case, where we solve the inverse problem by 
exploiting relevant results, purposely introduced through an auxiliary, direct formulation. In particular, 
we accurately report the solution method for the simply connected elliptical cross-section. Conversely, 
we omit details for the flow in the elliptical annulus between two confocal ellipses, since it can be derived 
from the previous case by a straightforward modification of the boundary conditions, as mentioned in 
Section |2.2.2| corresponding numerical results are nonetheless reported in Section |4.2[ for the sake of 
completeness. 
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3.1. Flow in the circle: the inverse problem. Let us preliminarily introduce the following 
T-periodic functions: 



(3.1) A(t) = ]T A m e^ 4 , f(t) = ]T /^""' w(x,t) = J2 ^m{x)e^\ <J m = ^ 

m£Z ragZ mSZ 

The aforementioned works by Sexl [5T] and Womersley [S9] addressed the direct problem, that is with 
assigned pressure gradient X(t). Corresponding solution can be constructed as an infinite series of Bessel 
functions Jo of order zero (Womersley considered a single Fourier mode, yet superposition can be invoked 
under the assumptions in [5J, since the corresponding sequence formally converges). As regards the inverse 
problem, an explicit map between the Fourier coefficients of the velocity w m and those of the flow rate 
f m has been recently obtained by the authors [5] . Relevant relations read 



irR 2 



Fi(;2 ;i Wo^ m /4)\ A m ^ _ / Jo((-l) 3/4 Wo„„)\ \ T 



oFi (;l;iWo^ m /4)y iw m ' m \ J ((-l) 3 / 4 Wo fl , m ) J iu> m ' 

where qFi(; ■;■) denotes the regularized confluent hyper-geometric function, and 

(3.2) Wo r>m = r s j^, 

is a non-dimensional parameter which generalizes the classical Womersley number. 

3.2. Flow in the ellipse: the auxiliary direct problem. The explicit solution for the direct 
problem in an elliptical vessel was originally derived by Khamrui 32 (although some details were missing) 
and Verma |56) . following the same approach of Womersley. In particular, given a single harmonic 
pressure gradient A = e lUmt , with m € Z, the solution (in elliptical coordinates) is the following: 

oo . 2 

(3.3) u(r),0) = y]C 2n Ge 2n ('n> -<?)ce 2 „(0, -?), with q = m . 

^— ' 4 v 

n=0 



The functions Cei n and cei n in (3.3 1 respectively denote the ordinary and modified Mathieu func- 
tions [40], while C2n represent suitable constants, determined from the no-slip (i.e. Dirichlet) boundary 
condition. It must be remarked that, despite the aforementioned functions were originally introduced by 
Mathieu in the 19th century (to study the vibration of an elliptical membrane, and they are still applied 
in the field |12j). they still deserve attention from a computational viewpoint, since their evaluation is 
prone to numerical instabilities [S3]. Moreover, in order to get a ^-periodic solution through the Mathieu 
functions (as in the present case), one needs to evaluate the eigenvalues of the Laplacian: this involves 
infinite tridiagonal matrices and it is very expensive. Furthermore, practical limitations on the cllipticity 
e — /3/a arise when directly using these special functions; see e.g. [26) . where numerical instabilities are 
reported for e < 0.3 (besides interesting expressions for e — >• 1, that is for elliptical cross-sections degen- 
erating to the circular geometry) . In light of these points, the quest for robust numerical methods able 
to efficiently compute the Mathieu functions is still an open issue, and positive contributions are being 
provided. In particular, a stepwise procedure is proposed in [24], suitable for evaluating the Mathieu 
functions in correspondence to large complex arguments typically associated with viscous flows. Such a 
strategy exploits a proper blend of backward and forward recurrence techniques in order to enhance the 
convergence of the involved computations. 

In order to circumvent the aforementioned limitations, we propose a novel numerical strategy, in the 
spirit of Fourier analysis, which only involves the Mathieu functions in an indirect way. In particular, we 
extend the approach introduced in [35] for the stationary case to the time-dependent one: our approach 
is based on a Fourier analysis in the variables 9 and t, while the problem in the variable 77 is kept in the 



physical space. More in detail, by applying the change of variables introduced in (2.5), we firstly recast 
the time-dependent problem as follows: 

(3.4) u t (t, v , 0) J(,, 0) - „ + = Jir,, 0) X(t). 



Then, since we look for T-periodic solutions (besides needing 27r-periodicity with respect to 9), we make 
the following ansatz of separation of variables for the velocity: 



(3.5) 



u(t,r),6) 



m,n£Z 



inO iujmt 



hereafter considered in place of the corresponding one in (3.1 ). Furthermore, since we consider real- valued 

solutions, the following conditions immediately apply: A_ m = A m , /_ m = f m and U- m ,-n — "777,71 (latter 
relation implies, in particular, $s(u(rf)o,o) — 0). 

We then start by deliberately addressing a direct problem, hence by assuming a given pressure 
gradient. In order to stress the fact that we temporarily turn our attention to a direct problem, we 
introduce a minor notation change (just within this section) by denoting the pressure gradient as a(t) 

(3.6) 



. By plugging (3.5) into (3.4), we arrive at the following problem: Solve, for each meZ, 



n£Z 



iLJ m u m>n (r]) e inf) J(r),9) - v{u" mn {r]) - n 2 u m , n (r])) e in8 = J{rj,6) o m 

u m , n {b) = n£Z, 



with (r/,9) £ E' and E' defined as in (2.7). The considered problem can be immediately simplified by 
invoking symmetry: velocity must be unchanged by the transformations 9 t— > —9 and 9 t— > 7r — 9, i.e. by 
reflection over the ellipse axes of symmetry. This, by also recalling the above conditions on conjugacy, 
provides 



(3.7) 



for n odd and u r 



for n even, m £ Z. 



Hence, we only need to solve the equations for the Fourier modes u m ,n, with n,m £ Z both greater or 
equal to zero, with n even. This permits to reduce by a factor eight the computational burden. That 



said, by plugging into (3.6) the Fourier expansion of J given in (2.6) and by equating the corresponding 



even n-modes (at fixed m), we obtain the following infinite family of ordinary differential equations 
(3.8) 
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m,2n+2 
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where we use u m , n in place of u m ^ n (rj) for brev ity, and Wo ajTO is a generalized Womersley number based 
on the semi-focal distance a, defined as in (3.2 1. Please also notice that q = 1W0 2 m /4, with q introduced 



in (3.3) (cf. e.g. O [26]). Furthermore, as in the stationary case, we need to assign an additional 
condition at r\ = (where the change of coordinates becomes singular) in order to have a well-posed 
mixed boundary value problems for u m n . In particular, for the solution to be smooth also near r\ — 0, it 
is suitable to enforce symmetry for the velocity-profile, i.e. a vanishing velocity gradient normal to the 
inter-focal line. This reads (cf. [13 Eq. (9)]) d v u(0,6,t) = 0, that is «'(0) mi „ = 0. 



The equations (3.8) represent an infinite system of non-homogeneous Mathieu ordinary differential 
equations. In order to get a computable system, one main idea is to approximate such a system with 
a finite dimensional (large enough) system of coupled ordinary differential equations. However, this 
cannot be achieved by simply neglecting all equations which involve modes u m .n for large enough |n|, 
since there is an infinite coupling (indeed, for each n £ Z, the equation for n-th mode is coupled with 
those of the two closest modulo 2 n-modes, thus providing a tridiagonal system) and such a rough 
simplification would lead to a system without solutions. The main idea in the present study is that 
higher frequencies are asymptotically small, since they are terms of a convergent Fourier series. The 



fact that the series converges derives from the abstract results in [5J [ST] , which can be even simplified 
in the case of an assigned smooth and periodic pressure gradient. Hence, by assuming that f m vanishes 
fast enough as \m\ — > +00 (and this is more than reasonable for realistic flows, see also Section [4]), the 
assumed smallness of the higher modes is fully justified. Therefore, once chosen a cut-off index TV e N 
(N > 2 to keep at least the first equation having zero as right-hand-side), we neglect u m . n for all neZ 
such that |n| > 2N since, in light of the aforementioned assumption, this does not affect the solution in 
a significant manner. From a practical viewpoint and based on the symmetries highlighted above, this 
implies to drop off u m ^N+2 in the differential equation satisfied by u m! 2Ar, and to only solve the equations 
for 

u m,n with n — 0, 2, ... , 2/V . In this way, for each m E N U {0}, we end up with a closed tridiagonal 
system of linear ordinary differential equations: it naturally represents the differential counterpart of 
the tridiagonal matrix which is obtained when evaluating the eigenvalues of the Mathieu functions by 
truncating the corresponding algebraic linear system (see e.g. [T| I24j). 

Before stating the tridiagonal system we actually address, we remark that, by virtue of linearity, it 
is convenient for our purposes to directly assume a m = 1. Moreover, we denote the unknown by v m ^ n 
in order to stress that, besides adopting a unitary pressure gradient, we are approximating u myTl by 
truncation. In light of the above arguments and derivations, for each m E N U {0} we have to solve the 
following boundary value problem: 
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It is worth remarking how, instead of truncating an infinite dimensional algebraic linear system in 
order to evaluate corresponding eigenvalues (as in [24]), we performed a Galerkin-type truncation of 
the system of ordinary differential equations which stems from the Fourier separation of variables. It 
is in this sense that we only indirectly referred to Mathieu functions, thus circumventing the numerical 
challenges associated therewith (indeed, we end up with a very simple tridiagonal system, which can be 
straightforwardly solved). 

Remark 3.10. An alternative pathway for solving the considered direct problem would be that of \2V$ : 
By considering real (i.e. sine and cosine) Fourier series functions, to solve the following system for the 
unknowns (4> m ,'>pm)- —n-ipm — v^4>m + lj n<f) m — i/Aip m , coupled with vanishing Dirichlet conditions. 
Indeed, the above system can be used to prove the existence of time-periodic solutions in an abstract way. 
However, in view of our main objective of obtaining a benchmark solution for the flow in elliptical cross- 
sections, it seems less convenient to follow such a strategy, since it involves the bi-Laplacian operator, 
which seems not to be easily manageable by separation of variables in elliptical coordinates. 



In a practical framework, it is necessary to only solve the system (3.9 1 for m <G A4*, where M* = 



{0, 1, . . • , M*}, M* denoting a sort of upper bound for the harmonic content of the given flow rate (i.e. 
its highest yet still relevant t-frequency is assumed to be u>m*)- Moreover, once fixed the fluid (i.e. v) 



and the size of the cross-section (i.e. a and b), the solution of (3.9) is only affected by the adopted cut-off 
index N. Let us stress this point by introducing the more descriptive notation Vm%\ as an alternative 
Hence, it is mandatory to choose N large enough to get a proper approximation for all the 



to V m ,r, 

modes 



■J2N) 



in particular for m up to M*. Let us denote by N* the sought value of N, to remind 
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that it also depends on M*. Two assets can be introduced to the purpose: accuracy (magnitude of the 
cut modes must be negligible compared to that of the kept ones) and independence from N (truncation 
must only affect the kept modes in a negligible way). As discussed above, higher frequencies are expected 
to be small than lower ones (by Fourier convergence); hence the accuracy asset can be formulated by 
introducing the following non-dimensional quantity: 



(3.11) 



(M(m,N) 



J m, 2AM|oo 
v m,0 II oo 



where H^Hoo = max^g^v |u(?7)|, E' — [0,6] denotes the considered 77-range, and the n — mode is used 
to get a reference value. (Such a reference was used for the numerical tests in Section |4j turning out 
to be non- vanishing; otherwise, alternative choices can be introduced by exploiting larger n values.) In 
particular, once chosen a threshold ft, let us define an index TV? as the smallest integer N such that 
fj,(m,N) < p,, for all m E M* . Hence, by choosing N > we are guaranteed that truncation only 
cuts negligible n-modes, for all relevant values of m. As regards independence from N, the following 
non-dimensional quantity is introduced: 



(3.12) 



s(m, N) = max 



\v (2N) 



-{2N+2) 



\v {2N) \\ 

\ u m,0 Woo 



where Af = {0, 2, . . . , 2N — 2} (index N is not considered, for obvious reasons). Then, once chosen a 
threshold s, let us define an index N± as the smallest integer N such that s(m, N) < s, for all m E M*. 
Hence, by choosing N > Ng we are guaranteed that truncation only negligibly affects the kept n-modes, 
for all relevant values of m. Finally, we naturally define N* — max (N±, N± ) , so as to achieve both 
targeted assets. 

3.3. Flow in the ellipse: the inverse problem. The interest in a benchmark solution for the 
inverse problem in elliptical cross-sections has been highlighted in Section]!] Nevertheless, some relevant 
results have been only recently provided in |24j , grounded on the solution of the direct problem by direct 
computation of the Mathieu functions, as discussed in Section |3.2| (earlier attempts in [33] seem to be 
loosely grounded on a theoretical basis). In order to circumvent the aforementioned challenges related 
to such a computation, we propose an alternative method which is based on the fast-Fourier approach 



introduced in Section 3.2 Thanks to linearity, we simply need to link the Fourier coefficients of the flow 



rate with those of the pressure gradient (as in the stationary case, thus generalizing (2.8)). The relevant 
solution strategy is described below. 



The flow rate associated with the ansatz (3.5) reads 



mez 



m ez 



/ «m,nO?) ( / 
neZ Ja \Jo 



e ind J(r), 9) d9 drj 



Then, once substituted J from (2.6), by explicit calculations we get 

-tt/2 



1 



AnB 



J{r), 9)d9= < 



n = ±2, 

7rcosh(277) n = 0, 

neZ\{0,±2}, 
so that the following relation can be easily obtained, for each meZ: 



(3.13) 



fm 



2 rb 

aw 



- u m -2(r)) + 2cosh(2?7)M mi o(n) - w m , 2 (n) dr] 



Clearly, without an explicit knowledge of u m ,n(v) (i- e - °f the sought solution), it is not possible to 
evaluate the flux by means of (3.13). Nevertheless, by recalling from Section 3.2 the approximate 



solution v m ^ n (associated with a unitary pressure gradient), it is straightforward to exploit (3.13) in 
order to approximate the Fourier coefficients A m (of the unknown pressure gradient) in terms of f m . 
Indeed, thanks to linearity and by also recalling ( |3.7[ ), the following relation is immediately obtained, for 
each m € NU{0}: 

-1 



(3.14) 



A, 



fn 



(cosh(2n) Vmfiiv) - v m ,i{rj)) d V 
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Latter expression deserves a remark: the bracketed quantity must be non- vanishing for such an expression 
to be meaningful. If we consider the exact solution u m ^ n of the infinite dimensional system (based on 
explicit knowledge of the asymptotic behavior of appropriate Hilbcrtian norms of eigenfunctions of the 
bi-Laplace operator in the domain under consideration, as for the existence results proved in [5ll21j). then 
the map between f m and A m is one-to-one. This, in turn, implies a non-vanishing denominator in the 



expression corresponding to (3.14). Yet such a condition is not perfectly guaranteed when considering 
an approximate solution v m _ n . Nonetheless, if the numerical approximation is accurate enough (i.e. for 
a cut-off index N large enough), the approximate denominator is close enough to the exact one for the 
non-vanishing condition to be fulfilled. It must be pointed out that the same potentially critical issue 
occurs when directly using the Mathieu functions. For instance, in [24] the authors implicitly rely on 
a well-defined mapping between flow rate and pressure coefficients for obtaining their solution. In that 
case, the degree of the involved Mathieu functions must be large enough for the associated numerical 
method to be well-posed, in complete analogy with the present case. 

Having introduced the main ingredients, it is now possible to summarize the basic steps we propose 
for solving the inverse problem: 

(50) We start by choosing a fluid (thus the kinematic viscosity v) and the size (namely a and b) of 
the considered vessel cross-section. We also assume a given flow rate fit) having period T; 

(51) We determine an integer M such that the Fourier spectrum of f(t) is suitably approximated by 
considering M modes; popular metrics (e.g. the Pearson correlation coefficient used in (4U1 |2~4"] ) 
can be used to the purpose. Practical constraints, due for instance to the experimental procedure 
implemented for measuring /, may come into play at this stage (for instance, the Nyquist 
limit is considered in [24] ) . Anyway, by the end of this step the Fourier components f m , for 
m — 0, 1, . . . , M, are known; 

(52) We fix an index M* such that M* > M and, once fixed desired (small enough) thresholds p, and 



s, we determine TV* as detailed at the end of Section 3.2 By the end of this step, the (accurate 
and truncation-independent) computation of the modes u m .„ is achieved, for m = 0, 1, . . . , M* 
and n = 0, 2, . . . , 2N*; 

(S3) By considering the modes v m , n with m = 0, 1, . . . , M (and n — 0,2, ... , 2N*), we compute A m 



by means of (3.14|, by also exploiting the relevant coefficient f m determined at step (SI). By 



the end of this step, the coefficients X m are available, for m = 0, 1, . . . , M; 



(S4) We finally compute the sought approximate solution by replacing (3.5 1 with the following sum- 
mation: 

(3.15) v(t,n,B)Si X m cpm(v,0)e iu '- t , with Vm (tj, 9) = £ v m<n {<n) e m6 , 

m£M neJV"* 

where M = {—M, —M + 1, . . . ,M - 1, M} and Af* = {-2N*, -2N* + 2,..., 2N* - 2, 2N*}, 
having previously obtained the modes associated with negative values of m and n by conjugacy, 
as detailed in Section [3~2"1 

We conclude by adding two remarks on Step (S2). Firstly, to accomplish such a task, we need to solve 



the system (3.9) several times, namely for several N. However, this burden is not specifically introduced 
by our strategy: for instance, to iterate computations up to reaching an adequate level of accuracy is 
also compulsory when directly using the Mathieu functions (see e.g. [24]). Indeed, our strategy (which 
simply involves the solution of tridiagonal systems) holds potential for a more computationally-efficient 
procedure than direct management of the Mathieu functions. Latter remark regards the introduction of 
M*. In particular, if one is interested in a single flow rate it makes sense to directly choose M* = M , so 
as to minimize computations. However, it can be of interest to also address multiple flow rate conditions, 
e.g. to identify the effect of specific m-harmonics on the solution (this is of interest, for instance, for 
the development of pulsatility-based devices, as mentioned in Section [l]). In such a case, it is convenient 
to choose M* > M when taking step (S2), in particular by adopting an M* large enough to account 
for all the foreseen m-harmonic contributions. Indeed, once built the corresponding v m ^ n "basis" , it is 
possible to immediately explore many solutions, by simply iterating steps (S3)-(S4) in correspondence 
of any given set of f m coefficients. This leads to a very computationally-cheap procedure, since it is 
only required to evaluate some integrals, besides trivial summations. This aspect further supports the 
effective exploitability of the solution strategy we propose. 



4. Numerical results. The numerical strategy detailed in Section |3.3| is hereafter exemplified, 
based on measured data of blood and cerebrospinal fluid (CSF) flow under physiological conditions. 

11 



As anticipated in Section [T] despite the exploitation of physiological flow conditions, the considered 
test-cases do not lay strong claims of being physiologically representative, due to the assumption of a 
fully-developed flow, which can be hardly fulfilled in real situations. Main aim of such simulations is 
to quantitatively assess the gain in computational efficiency which is enabled by the proposed Fourier- 
based numerical solution, as compared to more elaborate finite element solvers. This, in turn, renders 
our benchmark solution effectively usable for developing more ambitious numerical approaches, up to 
serving as a debugging tool for complex 3D codes. 

We point out that, besides the simulations described below, which are associated with mild eccen- 
tricity values (above 0.5), we successfully manage to also apply our numerical strategy to more idealized 
geometries with eccentricity e.g. below 0.1, thus succeeding where previous approaches encountered dif- 
ficulties |26j . However, corresponding numerical results are very similar to those reported in the sequel 
and therefore they are not reported, for ease of presentation. 



4.1. Blood flow in the internal carotid artery. A flow rate waveform for blood flow within the 



human internal carotid artery (ICA) was adapted from [29J . see Fig 4.1 (left); corresponding period and 
period-averaged flow rate are T — 0.95 [s] and f = 4.11 [cm 3 /s], respectively. An average radius for such 
a vessel is 0.25 cm [55] . Hence, by assuming e.g. an eccentricity equal to 0.6, we consistently introduced 
an elliptical cross-section with semi-axes a — 0.25 [cm] and j3 = 0.15 [cm], being the corresponding 
77-range E' v = [0,6], with b — 0.69 (a more realistic vessel deformation, e.g. due to compression by the 
surrounding tissue, might imply some larger value for a; however such a refined estimate is beyond the 
scope of the present simulations). Furthermore, a characteristic diameter for the considered section is 
5 = a + (3 = 0.4 [cm], while the sect ion- averaged speeds is w = fa/ A = 34.9 [cm/s], where A = naf3 
represents the cross-sectional area. The adopted flux was then approximated by M = 15 modes; corre- 



sponding Fourier sum, also reported in Fig. 4.1 suitably replaces the data at hand (associated Pearson 
correlation coefficient differs from 1 by less than 10~ 3 ). Moreover, once introduced a characteristic speed 
w = A~ x max t / TS [ .i] f(t) — 58.36 [cm/s] and by assuming v = 3.5 ■ 10~ 2 [cm 2 /s] [50], it is possible to 
label the considered blood flow with a Reynolds number Re = wb/v = 670. In addition, once defined a 
characteristic frequency uj — ^X)m=o Ifml^rnj I ^m=o l/mlj) it is possible to also introduce a charac- 
teristic Womersley number, namely Wo = (8/2) ^ui/v = 1.5. Despite the scarcity of results on stability 
for pulsatile flows (which seem to be often contradictory with one another already for the circular cross- 
section [53]), the above figures suggest that the flow at hand is laminar. This was derived by considering 
available experimental thresholds for the circular cross-section [11] [25], possibly defined through empirical 
fittings as in [42j (to the purpose, it may be appropriate to also define relevant Reynolds and Womersley 
numbers in a slightly different way |llj . yet we verified that this does not affect the above statement on 
stability for the considered flow). That said, more detailed considerations on laminar stability cannot be 
drawn in the present scope: they certainly require more consolidated acquisitions, as expressly remarked 
e.g. in [IT]. Step (S2) of the procedure proposed in Section 3.3 was then addressed, by choosing M* = 50 
for the sake of illustration. More in detail, the contours shown in Fig. 4.1 were firstly derived for both 




Fig. 4.1. Blood flow within ICA. (Left) Normalized flow rate f/fo-' data, adapted from J29( , are plotted against 
Fourier reconstruction based on M = 15 modes. (Middle) Contour plots for /j,(m,N): e.g. it is necessary to choose 
2N = 34 for the relative magnitude of discarded modes to be below 10~ 12 , for all the considered m values. (Right) 
Contour plots for s(m,N): e.g. it is necessary to choose 2N = 28 to get a relative sensitivity to truncation below 10 — 12 , 
for all the considered m values. 
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jU(m, AT) and s(m,iV), respectively defined in (3.11) and (3.12), for m = 1, . . . ,M* and iV large enough. 
To the purpose, we coded the boundary value problem (3.9) within Matlab environment and we solved 
it by a shooting technique. As a result, by choosing e.g. fl = s = 10~ 12 , it is immediate to get AT? = 17 
and N~ — 14 from the relevant plots in Fig. 4.1 so as to finally determine N* — 17. Please also notice 
how, for any fixed m, ||w m) n||oo > ||^m,n+2||oo for n = 0, ...,2iV— 2, thus confirming the expected relative 
less influence of higher n-modes. Relevant modes v m n were thus obtained by solving (3.9) with N = N*; 
Fig 4.2 shows some solution components, including tpi defined in (3.15). 



. 0.04 
0.02 



. -0.02 
-0.04 



3 ~ 2 ' 
U -4- 




(m=l,n=2) 



. 0.04 
0.02 



■ -0.02 
-0.04 



ST 4 r 
1 2 



— -0.05 

5 -0.1 

a 

-0.15 
0.2 



n [cm] -0.2 -0.5 



Fig. 4.2. Blood flow within ICA. Real part (upper row) and imaginary part (lower row) of selected solution compo- 
nents, computed with 2N = 34: vi^ivfte 12 ® (left column), «i,4(»j)e l4e (middle column) and <pi(rj,6) (right column). Latter 
quantity, defined in (|3.15[l, accounts in particular for all the m = 1 modes. 



We then compared our numerical results with those achieved by means of the commercial finite 
element (FE) solver ADINA 8.8.1 (ADINA R&D Inc., MA, USA), available to the group. Such a solver 
uses a standard Galerkin formulation (stabilized by upwinding for higher Reynolds numbers while 
time-advancing is implemented by a first-order Euler backward method. As for FE simulations, a pipe 
with length I = 1605 was defined, in order to obtain a fully-developed flow in the central portion of the 
domain (such a condition was a posteriori checked). Furthermore, the flow rate shown in Fig |4.1| was 
imposed (by software coding) at the inlet cross-section, the no-slip (i.e. Dirichlet) boundary condition was 
enforced on vessel wall, and a reference pressure value was imposed at the outlet cross-section (such value 
is immaterial, due to the incompressible formulation). In addition, as regards simulation set-up, both 
space- and time-discretization were incrementally refined, up to obtain discretization-independent results. 
In particular, pipe domain was discretized by nearly 8.2 • 10 5 second-order accurate brick elements, and 7 
pulsation periods were simulated (time-periodicity was obtained after 4 periods). Finally, all simulations 
were run on a single core of a PC with Intel Core i7-960 3.20 GHz CPU and 24 GB RAM. Obtained results 
are shown in Fig 4.3 as expected, a very good agreement is achieved between the considered approaches, 
main difference residing in computational times. Indeed, once accomplished in roughly 3 minutes the 
preliminary step (S2) (see Section 3.3 ) so as to get truncation-independent results, the computation is 
finalized in a very few seconds. Conversely, runtime for the FE simulation was about 9 days; additional 
days were necessary to set-up the simulation (namely to get grid-independence, periodicity and fully- 
developed flow conditions), so that the overall computational time for the FE run can be estimated over 
20 days. Corresponding speed-up factor is close to 10 4 . 

4.2. Cerebrospinal fluid flow in the spine. A flow rate waveform for cerebrospinal fluid (CSF) 
in the upper cervical region of the human spine |31j ) was adapted from [24] , see Fig 4.4 corresponding 
period-averaged flow rate is /o = —0.11 [cm 3 /s] (negative sign here indicates that it is directed towards 
the lumbar region). Corresponding period was assumed equal to that one of the cardiac cycle, reported 



in Section 4.1 The considered cross-section can be approximated by the annulus between two confocal 
ellipses with a-i = 1.11 [cm], = 0.93 [cm] (from which a — 0.61 [cm]) and /3i = 0.43 [cm], being 
the corresponding ?y-range EL = [61,62]; with bi — 0.66 and b 2 = 1.21. Furthermore, a characteristic 
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Fig. 4.3. Blood flow within ICA. Normalized velocity profiles w/vj along the major (left) and minor (right) semi-axes 
of the ellipse, for selected non-dimensional times t/T. Results obtained by the proposed Fourier-based approach (solid 
curves) are compared with those achieved by a commercial FE solver (filled circles). 



thickness for the considered section is f = (0:2 — ct\ + 02 ~ Pi) /2 = 0.43 [cm], while the section-averaged 
speeds is w — fo/A = —0.47 [cm/s], where A represents the cross-sectional area. Such a flux was then 
approximated by M = 15 modes; corresponding Fourier sum, also reported in Fig. 4.4 suitably replaces 
the data at hand (associat ed P earson correlation coefficient differs from 1 by less than 10 -6 ). Moreover, 



by reasoning as in Section 4.1 with f in place of 6, and by assuming v = 10 2 [cm 2 /s] [37], it is possible 



to label the considered blood flow with a Reynolds number Re — wf/v = 64 and a W omersley number 

the above figures 



4.1 



Wo = (f /2)Wui /v = 5.5. In light of the relevant references introduced in Section 
suggest that the flow at hand is laminar. However, also in this case more consolidated results are needed 
for a stronger statement on stability (available data for the circular cross-section [54] [IT] better describe 
almost purely oscillating flows like the present one, yet not in an annular section); such a study is clearly 
beyond the present scope. By proceeding as for the ICA, we then determined N* by considering M* = 50 
m-modes: once chosen ft — s — 10~ 12 we got N* = 17 (associated contours are not shown for brevity, 
being similar to those in Fig 4.1 ). Relevant modes v m , n were thus obtained by solving (3.9 ) with N — N*; 



GSF flow rate within spinal cavity 




Fig. 4.4. CSF flow in the spine. (Left) Normalized flow rate f/fo: data, adapted from [Sfi, are plotted against 
Fourier reconstruction based on M = 15 modes. (Right) 3D Reconstruction of the CSF domain, as obtained from patient- 
specific magnetic resonance imaging (MRI). Cranial sub-arachnoid space \31j and upper cervical spinal region are shown; 
location of the considered annulus is shown in the inset. 
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Fig. 4.5. CSF flow in the spine. Real part (upper row) and imaginary part (lower row) of selected solution components, 
computed with 2N = 34: «i,2(??)e l2fi (left column), ?i I 4(»?)e l4e (middle column) and (pi(r],d) (right column). Latter 
quantity, defined in (|3.15[l, accounts in particular for all the m = 1 modes. 



Fig 4.5 shows some solution components, including ipi denned in (3.151 



As for the ICA test-case, we compared our numerical results with those achieved by the commercial 



solver ADINA (relevant details can be recalled from Section 4.1 1. For such FE simulations, a pipe with 
length £ = 100r was defined, in order to obtain a fully-developed flow in the central portion of the 
domain (such a condition was a posteriori checked) . Furthermore, the flow rate shown in Fig |4.4| was 
imposed (by software coding) at the inlet cross-section, the no-slip (i.e. Dirichlet) boundary condition 
was enforced on vessel walls, and a reference pressure value was imposed at the outlet cross-section (such 
value is immaterial, due to the incompressible formulation). In addition, as regards simulation set-up, 
both space- and time-discretization were incrementally refined, up to obtain discretization-independent 
results. In particular, pipe domain was discretized by nearly 3.8- 10 5 second-order accurate brick elements, 
and 7 pulsation periods were simulated (time-periodicity was obtained after 2 periods) . Obtained results 



are shown in Fig 4.6 (they are obviously identical to those reported in [24]). As expected, a very 
good agreement is achieved between the proposed Fourier-based solution and that provided by the finite 
element solver, main difference residing in computational times for the present case as well. Indeed, once 



accomplished in roughly 3 minutes the preliminary step (S2) (see Section 3.3) so as to get truncation- 
independent results, the computation is finalized in a very few seconds. Conversely, runtime for the FE 
simulation was about 4 days (less than for ICA thanks to the lower average speeds); additional days 
were necessary to set-up the simulation (namely to get grid-independence, periodicity and fully-developed 
flow conditions), so that the overall computational time for the FE run can be estimated over 10 days. 
Corresponding speed-up factor is well over 10 3 also in the present case. 

5. Concluding remarks. We proposed a novel numerical method for solving the inverse problem 
of pulsatile viscous flows in elliptical vessels and annuli. In particular, we addressed a fully-developed 
flow of an incompressible Newtonian fluid, so as preliminarily derive some analytical relations (otherwise 
hardly achievable) upon which to build our numerical solution. We are therefore aware of the fact that 
the solution we propose is hardly applicable to 3D finite-length complex domains (such a limitation, 
despite being obvious, is often understated by the many authors who also assumed a fully-developed 
flow). On regard, under quite general assumptions (yet still adopting an incompressible Newtonian fluid), 
many valuable approaches for problems with flow rate conditions have been proposed in the literature 
(see e.g. |18j). and they are the subject of ongoing research. Moreover, a fully numerical approach is 
the only viable when also considering non-Newtonian fluids, already when tackling stationary problems 
(possibly within curved or non-straight vessels as in [5]). Nevertheless, and in spite of the linear context 
within which it was derived, our solution to the considered inverse problem is original and non-trivial. 
Indeed, it provides an easily computable benchmark (as well as an approximation) for such a problem. In 
particular, by invoking separation of variables (in elliptical coordinates) we obtained a tridiagonal system 
of linear, ordinary differential equations and we obtained the sought solution by means of a Galerkin-type 
truncation. This method seem to theoretically corroborate the early attempts in |49] (in which proper 
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Fig. 4.6. CSF flow in the spine. Normalized velocity profiles w/vj along the major (left) and minor (right) semi-axes 
of the annulus, for selected non-dimensional times t/T. Results obtained by the proposed Fourier-based approach (solid 
curves) are compared with those achieved by a commercial FE solver (filled circles). 



references to the inverse problem are missing). Moreover, it represents an alternative to the numerical 
approach in [21], heavily based on a wise computation strategy of the Mathieu functions. Indeed, 
truncation of the aforementioned system represents the differential counterpart of the truncation of the 
infinite dimensional algebraic linear system in [24j , which is performed in order to evaluate corresponding 
eigenvalues. In this sense, we only indirectly referred to Mathieu functions, thus circumventing the 
numerical challenges associated therewith (see [23] HH1 [S3]). As a consequence, even if the considered 
differential system is not completely integrable (due to the coupling between Fourier modes which is 
implied by the elliptical cross-section already when addressing the stationary problem), our solution can 
be easily computed, even for small ellipticity values. Such an asset was shown through the numerical 
simulations reported in Section [4] based on measured flow rates for blood in the internal carotid artery 
(ICA), and cerebrospinal fluid (CSF) flow in the upper cervical region of the human spine. Indeed, despite 
being expected, the main advantage of the proposed method, as compared to fully numerical approaches, 
resides in computational efficiency: a speed-up over 10 3 was obtained for both the aforementioned test- 
cases. Such a gain is comparable with the one reported in |24j . even if simulation set-up times seem to 
be slightly understated therein (in particular, computation of the Mathieu functions was recognized as 
the most critical part, and several related convergence results are reported in such a valuable paper, yet 
at the same time corresponding computational burden was judged as inexpensive |24j ) . An additional, 
potential advantage of the proposed Fourier-based approach resides in the possibility of an extensive and 
detailed control over numerical implementation (as opposite to black-box libraries usually exploited for 
computing special functions), which could be further optimized e.g. by using consolidated techniques for 
Fast Fourier solvers [23] . 

These results effectively support the exploitation of the proposed method in a wider research scope, 
e.g. as a debugging tool and/or an improved source of boundary data for more ambitious numerical ap- 
proaches, possibly based on more realistic data. On regard and with reference to the study of blood flow, 
the proposed solution could provide an enhanced inflow condition for those approaches which merely 
scale a (steady) parabolic solution by a factor accounting for flow rate variability (see e.g. [55]). More- 
over, the proposed solution can be directly used for computing wall shear stress distribution along the 
adopted cross-section boundary. Indeed, vessel ellipticity affects such a stress [21], which in turn has been 
identified as a key player in the onset and development of relevant pathologies likes atherosclerosis, even 
if several acquisitions needs to be consolidated (see e.g. [3H] 10J- Furthermore, still within the context of 
bio-fluid dynamics, many investigations are being carried out for elucidating the role of CSF dynamics 
within human spine, in relation to the onset and progression of relevant diseases like syringomyelia. The 
proposed benchmark flow holds potential for application to such studies as well, in particular for assign- 
ing more accurate inflow/initial conditions in the cervical portion of the spinal fluid domain (which is 
usually described as a circular annulus between two finite-length, compliant cylinders [8] I14[ [16]: wave 
propagation is usually studied therein). It is worth mentioning that relevant investigations for such 
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fluid-structure interaction problems are being performed in the context of hemodynamics as well (see 
e.g. [3]). Furthermore, besides studies purely focusing on fluid dynamics, the proposed solution can be 
effectively applied to additional biomedical applications (as anticipated, the biomedical field is an elec- 
tive application domain, since flow rate is a main available physical datum), e.g. to targeted delivery 
problems. For instance, the effect of flow pulsatility on magnetic particle targeting was originally exem- 
plified in [B] , in order to highlight peculiar effects with respect to the stationary case [37] widely adopted 
when addressing particle targeting in blood capillaries. Such an enhancement straightforwardly applies 
to magnetoresponsive carriers at large, e.g. to the controlled navigation of MRI-guided microrobots in 
the vasculature |3J. Indeed, there is a growing interest in the robotic field for the development of mi- 
cro/miniaturized interventional tools to be navigated (e.g. through magnetic fields) within bodily fluids. 
These systems, currently defined as microrobots |41] despite their intrinsic simplicity (system compart- 
mentalization is indeed almost inapplicable at the micro scale) , hold potential for effective application to 
medical fields, and both blood flow and CSF flow were identified as elective carrying fluids. To further 
support this point, it is worth mentioning that, besides carrier transport along the mainstream, active 
microrobot navigation was also proposed, e.g. in the CSF flow. More precisely, magnetic navigation of 
helical microrobots was proposed for application in the spinal CSF [19] (being the helical shape inspired 
by the propulsion strategy of bacterial flagella, well suited for low- Reynolds- number regimes). Inter- 
estingly, relevant CSF velocity profile, together with corresponding fluidic actions on the microrobot, 
seem to be slightly overlooked in these studies (while attention is payed to other physical aspects, up 
to gravity, in spite of its minor effects at the micro scale [5S] )• This seems to be partly due to the lack 
of easily computable - yet physically derived - flow conditions; such an issue might be mitigated by our 
benchmark flow solution, thus contributing to better assess the feasibility of such long-term approaches 
in bodily fluids, even if still in a simplified context. Finally, an additional, remarkable application field 
has been recently identified, namely the one of energy harvesting from bodily fluids. Indeed, use of 
magnetic fields was proposed as a key powering strategy for most of the aforementioned microrobots, 
mainly in response to the impossibility of using on-board batteries at the micro scale (alternative ap- 
proaches based on harnessing the power of biological entities like bacteria were also proposed, yet their 
discussion is out of the present scope). In fact, powering (closely linked to actuation) appears as the 
real bottleneck currently hampering microrobotic applications, and energy harvesting holds promise to 
solve such a major issue. In this spirit yet in the wider context of implantable miniaturized electronic 
devices for biomedical applications, an implantable fuel cell has been recently proposed [45], powered 
through energy scavenging from the CSF. Relevant estimates are reported in the cited reference in order 
to assess powering effective sustainability (they address oxygen and glucose transport within the CSF 
domain, mainly in the brain region but also encompassing the spinal one). However, it is expressly 
recognized in [45 that a deeper insight into the actual CSF dynamics is required to accurately assess 
the performance of the proposed system. Such a quest further supports the development of modeling 
strategies for the CSF flow (yet such a point can be immediately extended to other fluid conditions), 
and the benchmark flow solution we propose, despite the inherent simplifications we repeatedly pointed 
out above, could contribute to take a leap ahead in many of the aforementioned applications areas. 

To conclude, while developing the proposed numerical method we deliberately sacrificed some degree 
of physiological representativeness, yet in the name of direct and efficient computability, and thus wider 
usability of the obtained benchmark flow solution. The many applications discussed above encourage to 
exploit our solution in an wide interdisciplinary context, so as to further push the corresponding scientific 
and technological frontiers. 
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